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We study the super-counter-fluid(SCF) states in the two-component hardcore Bose-Hubbard 
model on the square lattice, using the quantum Monte Carlo method based on the worm(directed 
loop) algorithm. Since the SCF state is a state of a pair-condensation characterized by (a'b) 7^ 
0, (a) — 0, and (b) — 0, where a and b are the order parameters of the two components, it is impor- 
tant to study behaviors of the pair-correlation function (aibla^bj). For this purpose, we propose a 
choice of the worm head for calculating the pair-correlation function. From this pair-correlation, we 
confirm the Kosterlitz-Thouless(KT) charactor of the SCF phase. The simulation efficiency is also 
improved in the SCF phase. 
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Since the pioneering work by Greiner et a/.[IJ, there 
have been experimental developments in the field of ul- 
tra cold atoms in optical lattices. On theoretical sides, 
the Bose-Hubbard model which is the effective model of 
such systems have been well studied. The exact ground 
state phase diagrams of the standard Bose-Hubbard 
model have been obtained by the quantum Monte Carlo 
method [||-[5|. Recently, multi-component systems have 
been realized experimentally For example, Catani 

et al. 0] trapped heteronuclear bosonic mixtures of 87 Rb- 
41 K in optical lattices. Due to the multi-degree of free- 
dom, several exotic phases have been predicted theo- 
retically in this system Q. One of such phases is the 
super-counter-fiuid(SCF) phase[l(|- This state is a pair- 
condensation where A-particles and B-holes(A and B rep- 
resent each component) are paired by the strong repulsive 
interspecies interaction. The SCF state is characterized 
by (a* 6) ^ 0, (a) = 0, and (b) = 0. In Ref [|, the ground 
state phase diagram of two-component hardcore Bose- 
Hubbard model at a commensurate filling is revealed by 
the quantum Monte Carlo simulations based on the worm 
algorithm [ll| . The phase diagram consists of a number of 
distinct phases and the SCF is one of them that appears 
in the strongly correlated region when the asymmetry 
between the two kinds of particles is weak. 

As mentioned above, the quantum Monte Carlo 
method is a powerful tool for investigating quan- 
tum many-body systems. The worm(directed loop) 
algorithm [ill LL2| is one of the most effective quantum 
Monte Carlo methods. In this algorithm, we can mea- 
sure off-diagonal quantities such as two-point correla- 
tion function (ajat). However, in the simulation of two- 
species particle systems, the pairing correlation func- 
tion (aib\ojbj) can not be measured by the conventional 
method. Instead, several kinds of stiffnesses have often 
been used to detect the SCF stated, 0|. In order to 
measure the pair-correlation function which is related to 
the order parameter more directly, we can use t he q uan- 
tum Monte Carlo method proposed in Refs[3, Il5j . In 
the method, we sample the canonical ensembles by up- 



dates based on a worm operator which is always local in 
the imaginary time. This leads to much better statis- 
tics for equal-time off-diagonal quantities including the 
pair-correlation function. Although this method needs 
unconventional update procedures, this is valid when we 
simulate the system of fixed particle numbers. Since we 
can work in the canonical ensemble, we do not need to 
adjust chemical potentials to obtain desired fillings. 

In this Letter, we present a simple technique of measur- 
ing the pair-correlation function by the quantum Monte 
Carlo method in the grand canonical ensemble. In this 
method, we can perform conventional update procedures 
which needs only slightly extension. By this method, 
we study properties of the novel SCF state. Since 
the SCF phase is characterized by the order parameter 
(a^b), we discuss behaviors of the pair-correlation func- 
tion {a,ib\aj)j). Furthermore, we demonstrate this kind 
of worm improves simulation efficiency. 

The model we considered here is the two-component 
hardcore Bose-Hubbard model on a squared lattice. The 
Hamiltonian is given by 

H = - ^ t a a}aj - ^ hf>\bj + ^Un al n hi , (1) 

where a\{a-%) and b\(bi) are the bosonic cre- 
ation(annihilation) operators on the site i for two 
species of boson(A and B), n a i = aja;, t & {%) is the hop- 
ping parameter of A(B)-bosons, and U(>0) represents 
the repulsive interspecies interactions. Each species are 
hardcore bosons, which means that two bosons of the 
same kind cannot occupy the same site. The summation 
is over the nearest-neighbor pairs and the system 
size is defined by N — L 2 . The periodic boundary con- 
dition is applied. In what follows, we consider the case 
of half-filling for each component, (n a i) = (nbi) = 1/2, 
where ( ) means the thermal average. 

We used the quantum Monte Carlo method based on 
the directed loop algorithm(DLA) [12, HB| • In this algo- 
rithm, d-dimensional quantum systems are mapped to 
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FIG. 1: (Color online) Logarithmic plots of C pa i r (r) (left panel) and C a (r) (right panel) at different temperatures. In the left 
panel, the dashed line represents r 1 ' 4 which is expected at the critical temperature in the KT transition. In the right panel, 
the inset is semi-logarithmic plots of the large r region where exponential decays can be observed. Error bars are drawn but 
most of them are smaller than the symbol size(here and the following figures). We take t a as units of temperature. 



(d+ l)-dimensional classical systems, using the Feynman 
path integral representation. The classical systems are 
called world-lines, because these (d + l)-dimensional sys- 
tems are composed of the d-dimensional classical systems 
and the imaginary time axis. In the world-line quantum 
Monte Carlo method, world-line configurations are sam- 
pled. 

The most characteristic feature of the DLA is update 
procedures. For updates, we create two discontinuous 
points called a worm in world-lines. In order to define 
weights of configurations which contain worms, we con- 
sider the Hamiltonian H — ?7aQa~?7bQb, where the source 
term Q a is given by 

Q a = J2 I dr{a\{T) + ai {T)}/2, (2) 

and Qb is defined likewise. In the above equation, r is 
the imaginary time, (3 represents the inverse tempera- 
ture, and ai(r) = e rH aie~ rH . -q a and -q^ are coefficients 
which can be chosen to optimize simulation efficiency. By 
these source terms, we can generate world-line configu- 
rations with multiple worms in both of the two sectors, 
A and B, at the same time. However, updating pro- 
cedures using these configurations are generally compli- 
cated. Instead, we often use configurations which has 
just one worm for updates. Let us call a worm that 
works on the A(B)-bosons an A{B)-worm. When we 
apply the conventional method to two-component sys- 
tems, one of A- worms or B- worms is created somewhere 
in the space-time stochastically. Then one of the dis- 
continuous points (head) moves stochastically along the 
world-lines, updating particle numbers. When the head 
meets the other discontinuous point (tail) again, the worm 
is annihilated. In this algorithm, we can compute the 
two-point correlation functions (a\aj) by counting the 



number of times the A-worm's head passes the position 
(r,r) = (n — rj,0) relative to the tail pj], [Hi, where Ti 
is the coordination vector of the site i. 

In order to measure the pair-correlation function 
(ai&Jajfrj), we introduce a new worm for updates. The 
worm is represented by the source term — r? pa irQpair, 
where Q pa j r is written by 

Qpair = E M4(tMt) + Oi(r)6|(r)}. (3) 

We call this kind of worm a pair-worm. In worm cre- 
ations, one of the three kinds of worms(A-worm, B-worm 
and pair-worm) is chosen with nearly equal probabilities, 
by setting ?] a = ?/b = ?7pair- The other procedures are 
the same as the conventional ones, including the way to 
measure the pair-correlation function. It should be em- 
phasized that creations of a conventional worm are still 
necessary to satisfies the ergodicity. We performed ex- 
tensive checks of the code against exact diagonalization 
on small lattices. By using the pair-worm, simulation 
efficiency in the SCF phase can also be improved as ex- 
plained in the latter part of this letter. 

To investigate properties of the SCF state, the equal- 
time correlation functions 

C a (r) = (a\a 3 ) (4) 
£7b(r) = (b\b,) (5) 
C pair (r) = (aibla]bj) (6) 

are observed. In Fig. [1] we plot these correlation func- 
tions as a function of distance r = \ri — rj\, at various 
temperatures in the case of t a — t-t,,2zt a /U = 0.6. In 
the square lattice, the coordination number is z = 4. To 
obtain the results of L — 64, each run contains 5 x 10 
Monte Carlo steps for measuring physical quantities. We 
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FIG. 2: (Color online) Temperature dependence of the corre- 
lation ratio Cpair(£/2)/C P air(£/4) for different system sizes. 
The inset is the finite-size scaling plot. 



performed 256 independent runs for estimating the sta- 
tistical errors. The t a = % case is simulated, because the 
SCF state stabilizes most and is easy to analyze. Since 
Cb(r) is trivially equal to C a (r) in this case, it is omit- 
ted here. At high temperatures, e.g. T/t a = 0.16,0.20, 
all correlation functions decays exponentially. When the 
temperature becomes lower (T/t a = 0.04,0.095), C pa i r (r) 
behaves as power law decay. On the other hand, C a [r) 
still decays exponentially. This is a clear evidence of the 
SCF phase where A-particles and B-holes cannot support 
supercurrent individually but must form pairs to do so. 

The power-law decay of the correlation function 
is the characteristic property of the KT phase (l9j. 
We can make it clearer by further analysis of 
Cpair(^)- Let's consider the correlation ratio defined by 
Cpair(£/2)/Cp a i r (£/4). This quantity is independent of 
the system size at and below the critical point like the 
Binder ratio. Fig. [5] shows the temperature dependence 
of the correlation ratio. It can be seen that the correla- 
tion ratio for different system sizes overlap at low tem- 
peratures. This is a sign of the KT transition. Then, 
we make a finite-size scaling analysis, using the scaling 
form for the KT transition C pa i r (L/2)/C pair (L/4:)(T) = 
/(L/exp(c/vtT - T KT )Aa))0, S3, where the critical 



temperature Tkt and the unknown value c is to be esti- 
mated. In the KT transition, the finite-size scaling analy- 
sis is generally difficult due to the the singular divergence 
and the logarithmic correction. However, the correlation 
ratio is known to be a good estimator for determining 
the critic al p oint, because the logarithmic correction is 
cancelled [20|. In the inset of Fig. [2] we plot the resulting 
finite-size scaling of the correlation ratio. From this anal- 
ysis, the estimates T K r/t a = 0.092(2) and c = 0.80(2) are 
obtained. 

By snapshots of the world-lines(Fig. [3]), we also obtain 
the physical picture of the SCF state as follows: Due to 
the strong repulsive interaction, each site is occupied by 
one boson(A-boson or B-boson) . This is a feature of Mott 
states. In order to gain the hopping energy, A-boson and 



FIG. 3: (Color online) Cross-sectional snapshot of world-line 
configurations in the SCF phase. The left and the right panels 
are world-line configurations for the A-bosons and B-bosons 
respectively. Parameters are t a = tb, U/t& = 15, /?t a = 20 
and N — 4 x 4. Solid and dashed lines stand for occupied and 
vacant sites respectively. The discontinuous points are the 
positions where bosons hop in the direction vertical to the 
paper. Arrows show one of the approximate pair-hoppings. 
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FIG. 4: (Color online) Comparison between the DLA without 
using pair-worms and the DLA which includes pair-worms by 
as a function of U. 



B-boson often hop to each other's site in bonds where an 
A-boson is occupied on one side and a B-boson on the 
other side. This pair-hopping is the origin of the pair- 
condensation. As a results, the net currents of A-bosons 
and B-bosons are in opposite directions. 

We also discuss simulation efficiency of our method. 
For this purpose, we show the comparison between the 
DLA which includes the pair-worm and the DLA which 
does not include the pair-worm in Fig. |U For the com- 
parison, the compressibility of A-bosons 



1 d(n a ) 
(n a ) 2 d[i a 



(7) 



is used, because this quantity showed the most remark- 
able difference of all measured quantities, as explained 
below. In the above equation, [i a means the chemical 
potential of A-bosons and K a is proportional to the flue- 
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tuation of the total particle number of A-bosons. In both 
simulations, we performed 1.8 x 10 4 Monte Carlo runs for 
both thermalization and measuring quantities. By the 
finite-size-scaling analysis based on the [d + 1) dimen- 
sional U{\) universality class [9], we obtained the quan- 
tum critical point (U/t a ) c = 10.8(2). As the interspecies 
interaction U is increased beyond this point, the system 
transits from two-species superfluid (2SF) phase to the 
SCF phase. We have observed that results by two meth- 
ods in Fig. 2 are consistent within error bars in the 2SF 
phase, but not in the SCF phase. 

To understand this discrepancy, the important fact is 
that the typical length of worm's movement correctly re- 
flects the correlation length in the DLApj}. In the SCF 
phase, updates by the conventional A(B)-worms are often 
very local, i. e. restricted in a very narrow region, because 
the correlation length of A(B)-bosons are short-ranged. 
On the other hand, the method with the pair-worm en- 
ables to update configurations globally even in the SCF 
phase, because the correlation length of a^b is divergent. 
When we measure the fluctuation of the total number of 
A(B)-bosons, global worm updates are necessary. This 
is because the total particle number can be changed only 
when the worm head goes across the periodic boundary 
of the imaginary axis and return to the worm tail. That 
is the reason why the compressibility K a incorrectly van- 
ishes in the method without the pair-worm and does not 
vanish in the method including the pair-worm. Other 
quantities such as stiffnesses make no difference between 
the two methods within error bars. However, the above 
result means that the method without the pair- worm may 



produce erroneous results in the SCF phase and needs 
longer runs. We must be careful when we are using the 
conventional method, because most quantities are close 
to correct equilibrium values as long as we do not pay 
attention to the compressibility of A(B)-bosons. 

In conclusion, we have proposed the simple QMC tech- 
nique for measuring the pair-correlation function. By this 
method, we have calculated several correlation functions 
including the pair-correlation function in order to confirm 
that the pair-condensation is characteristic and essential 
to the SCF phase. Furthermore, we demonstrated that 
the new worm also improves the simulation efficiency in 
the SCF phase. In this improvement of simulation effi- 
ciency, the important fact is that we have introduced the 
worm corresponding to the order parameter a^b whose 
correlation length is divergent at the criticality. Our re- 
sults suggest that suitable worms should be used accord- 
ing to what kind of correlation exits. This idea is ex- 
pected to be applied to general multi-degree-of-freedom 
systems where various type of order parameters exit. 
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